////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//  -----------------------        Ebrahim Foulaadvand, 01 June 2012   -------------------------          //
//                                                                                                        //
// The programme "AirDragAscend" evaluates the ascent time as well as the maximum height of an upward     //
// thrown object with velocity v0 (from the ground) in the presence of air force (quandratic in v) in     //
// the constant gravitational field g. It also computes the descent time and the hit velocity when the    //
// object returns to the ground. The Euler-Richardson algorithm is used.                                  //                                                                          //                                                                                                                 //
//                                                                                                        //                                                                                                         //
////////////////////////////////////////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>

using namespace std;

double  tau=0.01,g=9.8,m=1.,c2=0.01,v0,velocity,vhit=0.,ac=0.,tmid=0.,ymid=0.,
        vmid=0.,Tascent,TascentAN,Tdescent,TdescentAN,Ymax,YmaxAN,Y,delv0=1.;

// The variables are explained in the "FreeFall" prpogramme.

int N=100000,t,S=200.;

double a (double &v); // acceleration function.

main()

{

vector<double> v(N+1,0),y(N+1,0);

//arrays "y" and "v" store the object position and velocity at each time step.

ofstream fileTdiff("Tdiff c2=0.01.plt"); //output file which sketches the computed |Tascent-Tdescent| vs v0.
ofstream fileTdiffAN("TdiffAN c2=0.01.plt"); //output file which shows the analytic |Tascent-Tdescent| vs v0.
ofstream fileYmax("Ymax c2=0.01.plt"); //output file for the computed maximum height vs v0.
ofstream fileYmaxAN("YmaxAN c2=0.01.plt"); //output file for the analytic maximum height vs v0.
ofstream fileYmaxFree("YmaxFree.plt"); //output file for the maximum height vs v0 without air drag.


for(int s=1; s<=S; s++){ // loop over values of thrown velocity v0

v0=s*delv0;


// ----------- upward motion --------------------

y[0]=0.;  v[0]=v0; // initial condition

//------------------------ Euler-Richardson Method------------------------

velocity=v0; t=0;

while (velocity>0.01){

tmid=t*tau + 0.5*tau;
ymid=y[t] + 0.5*v[t]*tau;
vmid=v[t] + 0.5*tau*a(v[t]);

v[t+1]=v[t] + tau*a(vmid);
y[t+1]=y[t] + tau*vmid;

velocity=v[t+1];

t+=1;

}//end while

// "t" shows the number of timesteps elapsed until the object velocity becomes zero (within precision of 0.01).

Ymax=y[t]; // Ymax is the maximum height reached by the object in its upward motion.

Tascent=tau*(t+1); // Tascent is the computed ascent time.

TascentAN=sqrt(m/(c2*g))*atan(sqrt(c2/(m*g))*v0); // TascentAN is the analytic ascent time.


//cout<<"Tascent= "<<Tascent<<endl;getch();


//----------- downward motion ----------------


y[0]=Ymax;  v[0]=0.; // initial condition.

Y=Ymax; t=0;

while (Y>0.01){

tmid=t*tau + 0.5*tau;
ymid=y[t] + 0.5*v[t]*tau;
vmid=v[t] + 0.5*tau*a(v[t]);

v[t+1]=v[t] + tau*a(vmid);
y[t+1]=y[t] + tau*vmid;

Y=y[t+1];

t+=1;

}//end while

// "t" shows the number of timesteps elpased for the downward motion i.e. from the maximum height
// until the object hits the ground.

vhit=v[t]; // vhit is the computed hit velocity.

Tdescent=tau*(t+1); // Tdescent is the computed descent time.

double teta=sqrt(c2/(m*g))*v0;

double expr=1/(cos(atan(teta)));

TdescentAN=sqrt(m/(c2*g))*log(expr+sqrt(expr*expr-1)); // TdescentAN is the analytic descent time.

//cout<<"Tdescent= "<<Tdescent<<endl;getch();

teta=atan(sqrt(c2/(m*g))*v0);

YmaxAN=-(m/c2)*log(cos(teta));

fileTdiff<<v0<<"  "<<(Tdescent-Tascent)<<endl;

fileTdiffAN<<v0<<"  "<<(TdescentAN-TascentAN)<<endl;

fileYmax<<v0<<"  "<<Ymax<<endl;

fileYmaxAN<<v0<<"  "<<YmaxAN<<endl;

fileYmaxFree<<v0<<"  "<<(v0*v0)/(2*g)<<endl;


}// end s


cout<<"end."<<endl; getch();

return 0;

}

//----------------------------------------------------------------

double a (double &v) {


return ac=-g-(c2/m)*fabs(v)*v;

}